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ABSTRACT 

We develop a formalism for the identification and accurate estimation of the strength of struc- 
ture formation shocks during cosmological smoothed particle hydrodynamics simulations. 
Shocks not only play a decisive role for the thermalization of gas in virialising structures but 
also for the acceleration of relativistic cosmic rays (CRs) through diffusive shock accelera- 
tion. Our formalism is applicable both to ordinary non-relativistic thermal gas, and to plasmas 
composed of CRs and thermal gas. To this end, we derive an analytical solution to the one- 
dimensional Riemann shock tube problem for a composite plasma of CRs and thermal gas. We 
apply our methods to study the properties of structure formation shocks in high-resolution hy- 
drodynamic simulations of the Lambda cold dark matter (ACDM) model. We find that most 
of the energy is dissipated in weak internal shocks with Mach numbers M ~ 2 which are 
predominantly central flow shocks or merger shock waves traversing halo centres. Collapsed 
cosmological structures are surrounded by external shocks with much higher Mach numbers 
up to A1 ~ 1000, but they play only a minor role in the energy balance of thermalization. 
This is because of the higher pre-shock gas densities within non-linear structures, and the sig- 
nificant increase of the mean shock speed as the characteristic halo mass grows with cosmic 
time. We show that after the epoch of cosmic reionisation the Mach number distribution is 
significantly modified by an efficient suppression of strong external shock waves due to the 
associated increase of the sound speed of the diffuse gas. Invoking a model for CR acceler- 
ation in shock waves, we find that the average strength of shock waves responsible for CR 
energy injection is higher than that for shocks that dominate the thermalization of the gas. 
This implies that the dynamical importance of shock-injected CRs is comparatively large in 
the low-density, peripheral halo regions, but is less important for the weaker flow shocks oc- 
curring in central high-density regions of haloes. When combined with radiative dissipation 
and star formation, our formalism can also be used to study CR injection by supernova shocks, 
or to construct models for shock-induced star formation in the interstellar medium. 

Key words: Shock waves - intergalactic medium - galaxies: clusters: general - cosmology: 
large-scale structure of universe - cosmic rays - methods: numerical 



1 INTRODUCTION 

1.1 Structure formation shock waves 

Cosmological shock waves form abundantly in the course of struc- 
ture formation, both due to infalling pristine cosmic plasma which 
accretes onto filaments, sheets and haloes, as well as due to su- 
personic flows associated with merging substructures (Quilis et al. 
1998; Miniati et al. 2000; Ryu et al. 2003; Gabici & Blasi 2003; 
Pavlidou & Fields 2005). Additionally, shock waves occur due to 
non-gravitational physics in the interstellar and intracluster media, 
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e.g. as a result of supernova explosions. Structure formation shock 
waves propagate through the cosmic tenuous plasma, which is com- 
pressed at the transition layer of the shock while a part of the kinetic 
energy of the incoming plasma is dissipated into internal energy 
of the post-shock gas. Because of the large coUisional mean free 
path, the energy transfer proceeds through collective electromag- 
netic viscosity which is provided by ubiquitous magnetic irregular- 
ities (Wentzel 1974; Kennel et al. 1985). 

Cosinologically, shocks are important in several respects. (1) 
Shock waves dissipate gravitational energy associated with hierar- 
chical clustering into thermal energy of the gas contained in dark 
matter haloes, thus supplying the intrahalo medium with entropy 
and thermal pressure support. Radiative cooling is then required 
to compress the gas further to densities that will allow star forma- 
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tion. (2) Shocks also occur around moderately overdense filaments, 
which leads to a heating of the intragalactic medium. Sheets and fil- 
aments are predicted to host a warm-hot intergalactic medium with 
temperatures in the range 10^ K < T < 10^ K whose evolution is 
primarily driven by shock heating from gravitational perturbations 
breaking on mildly nonlinear, non-equiUbrium structures (Hell- 
sten et al. 1998; Cen & Ostriker 1999; Dave et al. 2001; Furlan- 
etto & Loeb 2004; Kang et al. 2005). Thus, the shock-dissipated 
energy traces the large scale structure and contains information 
about its dynamical history. (3) Besides thermalization, collision- 
less shocks are also able to accelerate ions of the high-energy tail 
of the Maxwellian through diflFusive shock acceleration (DSA) (for 
reviews see Drury 1983; Blandford & Eichler 1987; Malkov & 
O'C Drury 2001). These energetic ions are reflected at magnetic ir- 
regularities through magnetic resonances between the gyro-motion 
and waves in the magnetised plasma and are able to gain energy 
in moving back and forth through the shock front. This accelera- 
tion process typically yields a cosmic ray (CR) population with a 
power-law distribution of the particle momenta. Nonlinear studies 
of DSA have shown that a considerable part of the kinetic energy 
flux passing through shocks can be channelled into non-thermal 
populations, up to about one-half of the initial kinetic energy of 
the shock (Berezhko et al. 1995; EUison et al. 1996; Malkov 1998, 
1999; Kang et al. 2002). Note that CRs have sufficient momentum 
not to resonate with the electromagnetic turbulence in the shock 
front itself They hence experience the shock as a discontinuity, 
i.e. the CR population is adiabatically compressed by the shock 
(e.g., Drury 1983). 

Indeed, CR electrons have been observed in the intra-cluster 
medium (ICM) of galaxy clusters through their diffuse synchrotron 
emission (Kim et al. 1989; Giovannini et al. 1993; Deiss et al. 
1997). In addition to these extended radio haloes which show a 
similar morphology compared to the thermal X-ray emission, there 
have been extended radio relics observed in the cluster periph- 
ery (e.g., Rottgering et al. 1997) which might well coincide with 
merger shock waves as proposed by EnBlin et al. (1998). Some 
clusters have also been reported to exhibit an excess of hard X- 
ray emission compared to the expected thermal bremsstrahlung of 
the hot ICM, most probably produced by inverse Compton up- 
scattering of cosmic microwave background photons by relativis- 
tic electrons (Fusco-Femiano et al. 1999; Sanders et al. 2005). It 
has been proposed that a fraction of the diffuse cosmological y- 
ray backgroimd radiation originates from the same processes (Loeb 
& Waxman 2000; Miniati 2002; Reimer et al. 2003; Berrington & 
Dermer 2003; Kuo et al. 2005). 

To date, there are two diflterent scenarios explaining these 
non-thermal emission processes. (1) Reacceleration processes of 
'mildly' relativistic electrons (y - 100 - 300) being injected 
over cosmological timescales into the ICM by sources like radio 
galaxies, supernova remnants, merger shocks, or galactic winds, 
which all can provide an efficient supply of highly-energetic CR 
electrons. Owing to their long lifetimes of a few times 10' years 
these 'mildly' relativistic electrons can accumulate within the ICM 
(Sarazin 2002), until they experience continuous in-situ accelera- 
tion either via shock acceleration or resonant pitch angle scattering 
on turbulent Alfven waves (Jaff'e 1977; Schlickeiser et al. 1987; 
Brunetti et al. 2001; Ohno et al. 2002; Brunetti et al. 2004). (2) In 
the ICM, the CR protons have lifetimes of the order of the Hubble 
time (Volk et al. 1996), which is long enough to difltuse away from 
the production site and to maintain a space-filUng distribution over 
the cluster volume. These CR protons can interact hadronically 
with the thermal ambient gas producing secondary electrons, neu- 



trinos, and y-rays in inelastic collisions throughout the cluster vol- 
ume, generating radio haloes through synchrotron emission (Den- 
nison 1980; Vestrand 1982; Blasi & Colafrancesco 1999; Dolag & 
Enl31in 2000; Pfrommer & EnBlin 2003, 2004a,b). Cosmological 
simulations support the possibility of a hadronic origin of cluster 
radio haloes (Miniati et al. 2001). 

1.2 Hydrodynamical simulations 

Hydrodynamical solvers of cosmological codes are generally clas- 
sified into two main categories: (1) Lagrangian methods like 
smoothed particle hydrodynamics (SPH) which discretise the mass 
of the fluid, and (2) Eulerian codes, which discretise the fluid vol- 
ume. SPH methods were first proposed by Gingold & Monaghan 
(1977) and Lucy (1977) and approximate continuous fluid quanti- 
ties by means of kernel interpolation over a set of tracer particles. 
Over the years, SPH techniques have been steadily improved and 
foimd widespread applications in cosmological problems (Evrard 
1988; Hemquist & Katz 1989; Navarro & White 1993; Springel & 
Hemquist 2002). 

In contrast, Eulerian methods discretise space and represent 
continuous fields on a mesh. Originally, Eulerian codes employed 
a mesh which is fixed in space (Cen & Ostriker 1993; Yepes et al. 
1995) or adaptively moving (Pen 1998), while more recently, adap- 
tive mesh refinement (AMR) algorithms have been developed for 
cosmological applications (Berger & Colella 1989; Bryan & Nor- 
man 1997; Norman & Bryan 1999; Abel et al. 2002; Kravtsov et al. 
2002; Refregier & Teyssier 2002), which can adapt to regions of 
interest in a flexible way. 

Grid-based techniques offer superior capabilities for captur- 
ing hydrodynamical shocks. In some algorithms, this can be done 
even without the aid of artificial viscosity, thanks to the use of Rie- 
mann solvers at the cell-level, so that a very low residual niuneri- 
cal viscosity is achieved. However, codes employing static meshes 
still lack the resolution and flexibility necessary to tackle structure 
formation problems in a hierarchically clustering universe, which 
is characterised by a very large dynamic range and a hierarchy of 
substructure at all stages of the evolution. For example, techniques 
based on a fixed mesh are seriously limited when one tries to study 
the formation of individual galaxies in a cosmological volume, sim- 
ply because the internal galactic structure such as disk and bulge 
components can then in general not be sufficiently well resolved. 
A new generation of AMR codes which begin to be applied in cos- 
mology may in principle resolve this problem. However, a number 
of grid-based problems remain even here, for example the dynam- 
ics is not Galilean-invariant, and there can be spurious advection 
and mixing errors, especially for large bulk velocities across the 
mesh. 

These problems can be avoided in SPH, which thanks to its 
Lagrangian nature and its accurate treatment of self-gravity is par- 
ticularly well suited for structure formation problems. SPH adap- 
tively and automatically increases the resolution in dense regions 
such as galactic haloes or centres of galaxy clusters, which are the 
regions of primary interest in cosmology. One drawback of SPH is 
the dependence on the artificial viscosity which has to deUver the 
necessary entropy injection in shocks. While the parametrization 
of the artificial viscosity can be motivated in analogy with the Rie- 
mann problem (Monaghan 1997), the shocks themselves are broad- 
ened over the SPH smoothing scale and not resolved as discon- 
tinuities, but post-shock quantities are calculated very accurately. 
However, to date it has not been possible to identify and measure 
the shock strengths instantaneously with an SPH simulation. 
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Being interested in dynamical implications of CRs on struc- 
ture formation and galaxy evolution, one faces not only the prob- 
lem of the interplay of gravity and hydrodynamics of a plasma 
composed of CRs and thermal particles but in addition radiative 
processes such as cooling and supernova feedback. To date, AMR 
codes have not yet matured to the point that they can address all 
these requirements throughout a cosmological volume, although 
there are recent efforts along these lines (e.g. Kang & Jones 2005; 
Jones & Kang 2005). It would therefore be ideal if SPH codes for 
structure formation could acquire the ability to detect shocks reli- 
ably during simulations. Previous work on shock detection in SPH 
simulations (Keshet et al. 2003) was restricted to a posteriori anal- 
ysis of two subsequent simulation time-slices, which can then be 
used to approximately detect a certain range of shocks as entropy 
jumps. 



1.3 Motivation and structure 

This article seeks to close this gap in order to allow studies of the 
following questions. (1) The cosmic evolution of shock strengths 
provides rich information about the thermal history of the baryonic 
component of the Universe: where and when is the gas heated to 
its present temperatures, and which shocks are mainly responsible 
for it? Does the missing baryonic component in the present-day 
universe reside in a warm-hot intergalactic medium? (2) CRs are 
accelerated at structure formation shocks through diffusive shock 
acceleration: what are the cosmological implications of such a 
CR component? (3) Shock waves are modified by nonlinear back- 
reaction of the accelerated CRs and their spatial diffusion into the 
pre-shock regime: does this change the cosmic thermal history or 
give rise to other effects? (4) Simulating realistic CR profiles within 
galaxy clusters can provide detailed predictions for the expected 
radio synchrotron and y-ray emission. What are the observational 
signatures of this radiation that is predicted to be observed with 
the upcoming new generation of y-ray instruments (imaging atmo- 
spheric Cerenkov telescopes and the GLAST' satellite) and radio 
telescopes (LOFAR^ and extended Very Large Array)? 

The purpose of this paper is to study the properties of struc- 
ture formation shock waves in cosmological simulations, allowing 
us to explore their role for the thermaUzation of the pristine plasma, 
as well as for the acceleration of relativistic CRs through DSA. 
In particular, we develop a framework for quantifying the impor- 
tance of CRs during cosmological structure formation, including 
an accounting of the effects of adiabatic compressions and rarefac- 
tions of CR populations, as well as of numerous non-adiabatic pro- 
cesses. Besides CR injection by structure formation shocks, the 
latter include CR shock injection of supernova remnants, in-situ 
re-acceleration of CRs, spatial diffusion of CRs, CR energy losses 
due to Coulomb interactions, Bremsstrahlung, and hadronic inter- 
actions with the background gas, and the associated -y-ray and ra- 
dio emission due to subsequent pion decay. A full description of 
these CR processes and their formulation for cosmological appli- 
cations is described in EnBlin et al. (2006), while the numerical im- 
plementation within the SPH formaUsm is given by Jubelgas et al. 
(2006). In this work we provide a crucial input for this modelling: 
a formalism for identifying and accurately estimating the strength 
of structure formation shocks on-the-fly during cosmological SPH 
simulations. 



' Gamma-ray Large Area Space Telescope, http://glast.gsfc.nasa.gov/ 
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The paper is structured as follows. The basic cosmic ray vari- 
ables are introduced in Section 2. The formaUsm for identifying 
and measuring the Mach number of shock waves instantaneously 
within an SPH simulation is described in Section 3 for a purely 
thermal gas as well as for a composite plasma of CRs and thermal 
gas. The numerical implementation of the algorithm is discussed in 
Section 4. In Section 5, we compare shock tube simulations to ana- 
lytic solutions of the Riemann problem which are presented in Ap- 
pendices A and B. Finally, in Section 6, we perform cosmological 
non-radiative simulations to study CR energy injection at shocks, 
and the influence of reionisation on the Mach number distribution. 
A summary in Section 7 concludes the paper. 



2 BASIC COSMIC RAY VARIABLES 

Since we only consider CR protons^, which are at least in our 
Galaxy the dominant CR species, it is convenient to introduce the 
dimensionless momentum p = /'p/(mp cught). CR electrons with 
y < 100 experience efficient Coulomb losses such that their energy 
density is significantly diminished compared to the CR energy den- 
sity (Sarazin 2002). The differential particle momentimi spectrum 
per volume element is assumed to be a single power-law above the 
minimum momentum q: 



f(p)- 



dN 
dpdV 



■ Cp-^eip-q). 



(1) 



6{x) denotes the Heaviside step function. Note that we use an effec- 
tive one-dimensional distribution function f{p) = 4np^p^\p). The 
CR population can hydrodynamically be described by an isotropic 
pressure component as long as the CRs are coupled to the ther- 
mal gas by small scale chaotic magnetic fields. The differential CR 
spectrum can vary spatially and temporally (although for brevity 
we suppress this in our notation) through the spatial dependence of 
the normalisation C = C{r, t) and the cutoff q = q{r, t). 

Adiabatic compression or expansion leaves the phase-space 
density of the CR population unchanged, leading to a momentum 
shift according \.o p p' = (p/poY^^ p for a change in gas density 
from po to p. Since this is fully reversible, it is useful to introduce 
the invariant cutoff and normalisation q^ and Co which describe the 
CR population via equation (1) if the inter-stellar medium (ISM) 
or ICM is adiabaticaUy compressed or expanded to the reference 
density po. The actual parameters are then given by 



Po 



1/3 



^0 and C(p) ■■ 



Co. 



(2) 



These adiabatically invariant variables are a suitable choice to be 
used in a Lagrangian description of the CR population. 
The CR number density is 



"CR 



-i: 



dpm = 



Cq' 



or- 1 ' 



(3) 



' a-particles carry a significant fraction of the total CR energy. Neverthe- 
less, the assumption of considering only CR protons is a reasonable ap- 
proximation, since the energy density of a-particles can be absorbed into 
the proton spectrum. A GeV energy or-particle can be approximated as an 
ensemble of four individual nucleons travelling together due to the rela- 
tively weak MeV nuclear binding energies compared to the kinetic energy 
of relativistic protons. 
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provided, that a > 1. The kinetic energy density of the CR popula- 
tion is 



C Mr, cl. 



SCR 



Jo a — I 



a -2 3 - a 



+«'-''(VrT^-i) 



(4) 



where Tp(p) = (y/l + p^- I) nip c^^^j is the kinetic energy of a pro- 
ton with momentum p, and h) denotes the incomplete Beta- 
function which is defined by SAa,b) = - 0*"'df. The 
integral of equation (4) is well-defined if we assume a > 2. The 
CR pressure is 

Pc. = ^Idpm/sp 



s 



■2 3- 



(5) 



where j3 = v/cngu = pi -\Jl + p- is the dimonsionless velocity of 
the CR particle. Note that for 2 < a < 3 the kinetic energy density 
and pressure of the CR populations are well defined for the limit 
q —> 0, although the total CR number density diverges. 

The adiabatic exponent of the CR population is defined by 



rcR = 



dlogPcR 



dlogp 



(6) 



while the derivative has to be taken at constant entropy 5 . Using 
equations (2) and (5), we obtain for the CR adiabatic exponent 



rcR = 



aPcR 9£ ^ 9fcR 
dC dp dq dp 



P 
Pqr 



a + 2 2 , 
— -3. -A.) 



l a-2 3-a \ 



(7) 



Note that in contrast to the usual adiabatic exponent, the CR adia- 
batic exponent is time dependent due to its dependence on the lower 
cutoff of the CR population, q. The ultra-relativistic limit {q 00) 
of the adiabatic exponent, where ycR ~* 4/3, can easily be ob- 
tained by using the integral representation of the incomplete Beta- 
function and applying a Taylor expansion to the integrand. In the 
non-relativistic limit « 1 and a > 3), the adiabatic exponent 
approaches yea 5/3. This can be seen by evaluating the CR 

pressure in this limit, Pck = C q"'" and applying the defini- 

tion of ycR in equation (6). Considering a composite of thermal and 
CR gas, it is appropriate to define an eff'ective adiabatic index by 



Teff = 



dlogCPtt + PcR) 



dlogp 



_ yth f th + ycR PcR 

^"111 + ^'CR 



(8) 



3 MACH NUMBERS WITHIN THE SPH FORMALISM 

The shock surface separates two regions: the upstream regime (pre- 
shock regime) defines the region in front of the shock whereas the 
downstream regime (post-shock regime) defines the wake of the 
shock wave. The shock front itself is the region in which the mean 
plasma velocity changes rapidly on small scales given by plasma 
physical processes. All calculations in this section are done in the 
rest frame of the shock which we assume to be non-relativistic. 
This assumption is justified in the case of cosmological structure 
formation shock waves for which typical shock velocities are of 
the order of 10' km s"\ 



Particles are impinging on the shock surface at a rate per unit 
shock surface, j, while conserving their mass: 



(9) 



Here V\ and V2 indicate the plasma velocities (relative to the shock's 
rest frame) in the upstream and downstream regime of the shock, 
respectively. The mass densities in the respective shock regime are 
denoted by pi and pi - Momentum conservation implies 



Pi +Pxv\ = P2 +P2W2 



(10) 



where P, denotes the pressure in the respective regime i £ {1,2}. 
The energy conservation law at the shock surface reads 

,,2 

(11) 



(ei + Pi)p7' + ^=is2 + P2)P2 + y • 



Si denotes the internal energy density in the regime i e {1, 2}. Com- 
bining solely these three equations without using any additional 
information about the equation of state, we arrive at the following 
system of two equations: 



,2»^2„2 ^ (P2-Pl)PlP2 
P2 -Pl 



P2 
Pl 



piMfc, - — 

2t> + P, + P2 
2s, + P| + P, 



(12) 

(13) 



Here we introduced the Mach number in the upstream regime. 
Ml = vi/ci, which is the plasma velocity in units of the local 
sound speed ci = s/yPi/pi.'* 

3.1 Poly tropic gas 

Non-radiative polytropic gas in the regime ; 6 { 1, 2) is characterised 
by its particular equation of state. 



Si = — —rPi or equivalently P, = Pol— 

r - 1 \po 



(14) 



where 7 denotes the adiabatic index. This allows us to derive the 
well-known Rankine-Hugoniot conditions which relate quantities 
from the upstream to the downstream regime solely as a function 
of Mi: 

(r+i)M? 



P2 ^ 

Pl (y-l)M] + 2' 

P2 ^ 2yM\-(y-X) 
Pl 



(15) 



(16) 



(17) 



7+1 

T2 _ [2yA1^-(y-l)] [(y-l)A1^ + 2] 

Ti " {y+\fM\ 

In cosmological simulations using a Lagrangian description 
of hydrodynamics such as SPH, it is infeasible to identify the rest 
frame of each shock and thus M.\ unambiguously, especially in the 
presence of multiple oblique structure formation shocks. As an ap- 
proximative solution, we rather propose the following procedure, 
which takes advantage of the entropy-conserving formulation of 
SPH (Springel & Hemquist 2002). For one particle, the instanta- 
neous injection rate of the entropic function due to shocks is com- 
puted, i.e. dA/d?, where A denotes the entropic function A{s) de- 
fined by 



P = A{s)p' 



(18) 



Note, that the symbol c (sometimes with subscript) denotes the sound 
velocity. 
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and s gives the specific entropy. Suppose furtlier that the shock is 
broadened to a scale of order the SPH smoothing length fi,h, where 
fh~2 denotes a factor which has to be calibrated against shock- 
tubes. We can roughly estimate the time it takes the particle to pass 
through the broadened shock front as Af = fi,h/ v, where one may 
approximate v with the pre-shock velocity Vi . Assuming that the 
present particle temperature is a good approximation for the pre- 
shock temperature, we can also replace v\ with M\C\.^ 

Based on these assumptions and using AA] ^ ArdAi/d?, one 
can estimate the jump of the entropic function the particle will re- 
ceive while passing through the shock: 



Ai + AA 



= 1 + 



fhh dAi 



MiCiAi dt 



A2 PilpiX 

T = "5" ~ =fA{Mi), 

M Pi \P2} 



where 

/a(A1i) = — ; 

y + 1 



(7- l)M\ + 2 



{y+l)M\ 



(19) 



(20) 



(21) 



using equations (15) and (16). Combining equations (19) and (20), 
we arrive at the final equation which is a function of Mach nimiber 
only: 



[/a(7V1i)-1]A1i = 



fhh dAi 
c\Ai df 



(22) 



The right-hand side can be estimated individually for each particle, 
and the left-hand side depends only on A1|. Determining the root 
of the equation hence allows one to estimate a Mach number for 
each particle. 



3.2 Composite of cosmic rays and tliermal gas 

In the presence of a gas composed of cosmic rays and thermal com- 
ponents, equations (9) to (13) are still applicable if one identifies 
the energy density e, and the pressure P, with the sum of the indi- 
vidual components in the regime / € {1, 2}, 



Pi = ^'cRi + Plhi- 



(23) 
(24) 



The sound speed of such a composite gas is c, = ^fy^sJPTfph where 
Teff.j is given by equation (8). Note that in contrast to the single- 
component fluid, for the general case there is no equivalent to the 
equation of state (equation 14) in terms of the total energy density 
because of the additivity of both pressure and energy density. 
For later convenience, we introduce the shock compression ratio 
and the thermal pressure ratio js. 



Xs = — and 

Pi 



Pth2 
-f Ihl 



(25) 



While taking the equation of state (equation 14) for the ther- 
mal gas component, we assume adiabatic compression of the CRs 
at the shock*". 



PcRiXs and ecR2 = scri^^s • 



(26) 



' Extensions of this approach that apply to particles within the SPH broad- 
ened shock surface will be described in Section 4.1. 
* Due to their much larger gyro-radii and high velocities, CR protons 
should not participate in the plasma processes of coUisionless shock waves. 



Here we assume a constant CR spectral index over the shock which 
holds only approximately owing to the weak dependence of the CR 
lower momentum cutoff q on the density (equation 2). 

For the composite of therinal and CR gas, it is convenient to 
define the effective entropic function A^g and its time derivative. 



(27) 
(28) 



Aeff = (P,h + ^'CR)P"''=^ 

dt dt ^ ' 

The expression for the time derivative of the effective adiabatic 
function uses the approximation of adiabatic compression of the 
CRs at the shock. Using the same assumptions like in the non- 
radiative case, we estimate the jump of the entropic function for 
the particle on passing through the shock made of composite gas: 

Combining equations (12), (13), (26), and (29), we arrive at 
the following system of equations, 

fl{Xs,ys) = Xs[P2(x„}\) - Pij 

X [P2(xs,ys)(x,pir"'-^'-'-''' - PiPr'^'f 

f2ixs,ys) = 2s2(x„y,) + Pi + P2(x„y,) 

- Xs [2ei + P| + P2{x„ >s)] = 0. (31) 

The effective adiabatic index in the post-shock regime is given by 

rcR^'cR2(Xs) + ythysPm 



yeS,2iXs,ys} = 



(32) 



P2(Xs,ys) 

Given all the quantities in the pre-shock regime, we can solve 
for the roots x^ and y^ of this system of two non-linear equations. 
This system of equations turns out to be nearly degenerate for plau- 
sible values of pre-shock quantities such that it might be convenient 
to apply the following coordinate transformation: 



(xs,y,) (x„Zs) with Zs 



ys - Xs 



(33) 



The Mach number Mi and the jump of internal specific energies 
can then be obtained by 



(P2-Pl)Xs , 

Ml = A —57 and 

Plcf(Xs - 1) 

"2 y, , P(a 

— = — where u = -— . 

Ml xs (rth - i)p 



(34) 
(35) 



4 NUMERICAL IMPLEMENTATION 
4.1 Poly tropic gas 

Applying the algorithm of inferring the shock strength within the 
SPH formalism in a straightforward manner will lead to system- 
atically underestimated values of the Mach number for SPH parti- 
cles which are located within the SPH broadened shock surface: the 
proposed algorithm of Section 3 assumes that the present particle 
quantities such as entropy, sound velocity, and smoothing length 
are good representations of the hydrodynamical state in the pre- 
shock regime, which is not longer the case for particles within 
the SPH broadened shock surface. To overcome this problem, we 
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define a decay time interval Afjec = rnin[/A/!/(A1ic), Af„ax]> dur- 
ing which the Mach number is set to the maximimi value that 
is estimated during the transition from the pre-shock regime to 
the shock surface. At this maximum, the corresponding particle 
quantities are good approximations of the hydrodynamical val- 
ues in the pre-shock regime. We thus have a finite temporal res- 
olution for detecting shocks, which is of order the transit time 
through the broadened shock front. Note that Af^ax is just intro- 
duced as a safeguard against too long decay times for very weak 
shocks. In the case of cosmological simulations, which are con- 
veniently carried out in a computational domain that is comov- 
ing with the cosmological expansion, we redefine the decay time 
A(lna)(iec = min[/f(a)^,/i/(Alic), A(lna)niax] where a denotes the 
cosmic scale factor and H{a) = a/a is the Hubble function. An ap- 
propriate choice for the safeguard parameter is A(ln a)max = 0.0025 . 

Secondly, there is no universal value fh which measures the 
SPH shock broadening accurately irrespective of the Mach number 
of the shock, especially in the regime of strong shocks. We there- 
fore use the original algorithm (with //, = 2) only for estimated 
Mach nvunbers with Mest < 3, while for stronger shocks, we apply 
an empirically determined formula (calibrated against shock-tubes) 
which corrects for the additional broadening of strong shocks and 
smoothly joins into the weak shock regime: 

Mc^ = {aMM'^ + CM exp-^-'/3-) ^^^^^ (3^,) 

where om = 0.09, bM = 1-34, and cm = 1-66. These nvunbers may 
depend on the viscosity scheme of the SPH implementation. 



4.2 Composite of cosmic rays and tliermal gas 

Our formalism of inferring the jump conditions for a composite of 
cosmic rays and thermal gas yields the density jump, = P2/P1, 
and the thermal pressure jump at the shock, = Pm/Pm (Sec- 
tion 3.2). As described in the previous section (Section 4.1), the 
values for the estimated jump conditions are systematically under- 
estimated in the regime of strong shocks (Mi ^ 5) implying an 
additional broadening of the shock surface. Thus, we proceed the 
same way as above: using the value of the density jump x^, we 
derive the Mach number of the shock through equation (34) and 
recalibrate it for strong shocks. In addition, we use the decay time 
Afjjc as before in the thermal case to obtain reliable Mach number 
estimates. The post-shock density is then obtained by multiplying 
the stored pre-shock density with the density jvunp Xs. 

In the case of a thermal pressure jump at the shock }\, we de- 
cided not to derive another empirical formula but rather exploit CR 
physics at non-relativistic shocks. Since the CR population is adi- 
abatically compressed at the shock in the limit of strong shocks, 
the total pressure jump is nearly solely determined by the jump of 
the thermal pressure in the post-shock regime, i.e. we can safely 
neglect the contribution of CRs to the pressure jump. This assump- 
tion is justified as long as the CR pressure is not dominated by sub- 
relativistic CRs of low energy which is on the other hand a very 
short lived population owing to Coulomb interactions in the ICM. 
Thus, the thermal post-shock pressure for Mi ^ 5 is estimated as 

P..^ ^^-^--^^-^ ., (37) 

rth + 1 

where Mi is obtained by equation (36), and Pi denotes the stored 
total pre-shock pressure. 



5 SHOCK TUBES 

To assess the reliability of our formalism and the validity of our 
numerical implementation, we perform a sequence of shock-tube 
simulations with Mach nvunbers ranging from M = 1.4 up to 
M = 100. We use a three-dimensional problem setup which is more 
demanding and more realistic than carrying out the computation in 
one dimension. Here and in the following, we drop the subscript 
'1' of the pre-shock Mach number for convenience. By compar- 
ing with known analytic solutions, we are able to demonstrate the 
validity of our implemented formalism. 

There exists an analytic solution of the Riemann shock-tube 
problem in the case of a fluid described by a polytropic equation 
of state, s = Pliy - 1) (cf. Appendix A). Unfortunately, a com- 
posite of thermal gas and CRs does not obey this relation. Thus, 
we derive an analytic solution to the Riemann shock-tube problem 
for the composite of CRs and thermal gas in Appendix B. This an- 
alytic solution assumes the CR adiabatic index (equation 6) to be 
constant over the shock-tube and neglects CR diffusion such that 
the problem remains analytically treatable. 



5.1 Polytropic tliermal gas 

We consider eight standard shock-tube tests (Sod 1978) which pro- 
vide a validation of both the code's solution to hydrodynamic prob- 
lems and our Mach number formalism. We consider first an ideal 
gas with 7 = 5/3, initially at rest. The left-half space {x < 250) 
is filled with gas at unit density, P2 = 1, and P2 = {y - 1) lO', 
while X > 250 is filled with low density gas pi = 0.2 at low 
pressure. The exact value of the low pressure gas has been cho- 
sen such that the resulting solutions yield the Mach numbers M = 
{1.4, 2, 3, 6, 10, 30, 60, 100) (cf. Appendix A). We set up the initial 
conditions in 3D using an irregular glass-like distribution of a total 
of 3 X 10"* particles of equal mass in hydrostatic equilibrium. They 
are contained in a periodic box which is longer in x-direction than 
in the other two dimensions, y and z. 

In the left-hand panel of Fig. 1, we show the result for the 
case of the Mach number A1 = 10 obtained with the GArK3ET-2 
code (Springel 2005; Springel et al. 2001) at time t = 0.5. Shown 
are the volume averaged hydrodynamical quantities {p{x)), {P{x)), 
{Vx(x)), and {M{x)) within bins with a spacing equal to the inter- 
particle separation of the denser medium and represented by solid 
black lines. One can clearly distinguish five regions of gas with dif- 
ferent hydrodynamical states. These regions are separated by the 
head and the tail of the leftwards propagating rarefaction wave, 
and the rightwards propagating contact discontinuity and the shock 
wave. The overall agreement with the analytic solution is good, 
while the discontinuities are resolved within 2-3 SPH smooth- 
ing lengths. Despite the shock broadening, the post-shock quanti- 
ties are calculated very accurately. Our formalism is clearly able to 
detect the shock and precisely measure its strength, i.e. the Mach 
number M. The pressure quantity drawn is not the hydrod3Tiam- 
ically acting pressure of the SPH dynamics but P = (y - l)pu, a 
product of two fields that are calculated each using SPH interpo- 
lation. Thus, the observed characteristic pressure blip at the con- 
tact discontinuity has no real analogue either in the averaged x- 
component of the velocity (v^ix)) or in the averaged Mach number 
{M(x)). The A-component of the velocity (l',(x)> shows tiny post- 
shock oscillations which might be damped with higher values of the 
artificial viscosity in the expense of a broader shock surface. The 
leftwards propagating rarefaction wave seems to exhibit a slightly 
faster signal velocity compared to the sound velocity. This might be 
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Polytropic thermal gas: Composite of CRs and thermal gas: 
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Figure 1. Shock-tube test canied out in a periodic tliree-dimensional box wliicli is longer in ,v-direction tlian in the other two dimensions where a shock with 
the Mach number A1 = 10 develops. The numerical result of the volume averaged hydrodynamical quantities (p(.v)), (P{x)), {Vt(x)), and (Ai(x)) within bins 
with a spacing equal to the interparticle separation of the denser medium is shown in black and compared with the analytic result in colour Left-hand panels: 
Shock-tubes are filled with pure thermal gas (y = 5/3). Right-hand panels: Shock-tubes are filled with a composite of cosmic rays and thermal gas. Initially, 
the relative CR pressure is Xcr = PcR/Pth = 2 in the left-half space (jc < 250), while we assume pressure equilibrium between CRs and thermal gas for 
X > 250. 



attributed to the SPH averaging process which obtains additional 
information on the SPH smoothing scale. 

In the left-hand panel of Fig. 2, we show the Mach number dis- 
tributions weighted by the change of dissipated energy per time in- 
terval, (di(,i,/(dlog A1)) for our eight shock-tubes. The sharp peaks 
of these distributions around their expected values log M are appar- 
ent. This demonstrates the reliability of our formalism to precisely 
measure shock strengths instantaneously during SPH simulations. 
The bottom panel shows their integral, i.e. the change of dissipated 
energy per time interval, (wth). The rising dissipated energy with 
growing Mach number reflects the larger amount of available ki- 
netic energy for dissipation. 

We additionally calculate the shock-injected CR energy using 



our formalism of diffusive shock acceleration described in EnBlin 
et al. (2006). However, the injected CR energy ^ucRjnj) was only 
monitored and not dynamically tracked. For comparison, we also 
show the theoretically expected injected CR energy = 
^inj ("th)i where ^i„j is the energy efficiency due to difl"usive shock 
acceleration (cf. EnBlin et al. 2006, for details). The good compar- 
ison of the simulated and theoretically expected shock-injected CR 
energy demonstrates that our formalism is reliably able to describe 
the on-the-fly acceleration of CRs during the simulation. 
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Poly tropic thermal gas: 
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Figure 2. Mach number distributions weigiited by the change of dissipated energy per time interval, (dHth/(d log M)) for our eight three-dimensional shock- 
tubes. Left-hand panels: Shock-tubes are filled with pure thermal gas (y = 5/3). Right-hand panels: Shock-tubes are filled with a composite of cosmic rays and 
thermal gas. Initially, the relative CR pressure is Xcr = Pcr/^'ui = 2 in the left-half space, while we assume pressure equilibrium between CRs and thermal 
gas. Bottom panels: Shown are the change of dissipated energy per time interval, <Mth> (shown with x), the shock-injected CR energy <McR,mj> (+). and the 
theoretically expected injected CR energy <m*^°^> (o) which is calculated following EnBUn et al. (2006). 



5.2 Composite of cosmic rays and thermal gas 

Again, we consider eight shock-tube simulations containing a com- 
posite of cosmic rays and thermal gas, providing a useful validation 
of our CR implementation in solving basic hydrodynamic problems 
as well as our Mach number formalism in the presence of CRs. 
In these simulations, we neither inject shock-accelerated CRs nor 
consider CR diffusion: these processes would lead to CR modified 
shock structures and shall be the subject of a companion paper. 

To characterise this composite fluid, we define the relative CR 
pressure Xqk = Pcr/Pa- Our composite gas is initially at rest, 
while the left-half space (x < 250) is filled with gas at unit den- 
sity, p2 = 1, XcR2 = 2, and Pthi = (y - I) 10^ while x > 250 
is filled with low density gas pi = 0.2, Xcri = 1, at low pres- 
sure. The exact value of the low pressure gas has again been 
chosen such that the resulting solutions yield the Mach numbers 
At = (1.4, 2, 3, 6, 10, 30, 60, 100) (cf Appendix B). Otherwise, we 
use the same initial setup as in Section 5.1. This CR load represents 
a rather extreme case and can be taken as the limiting case for our 
Mach number formalism in the presence of CRs. Cosmologically, 
it may find application in galaxy mergers where the outer regions 
might be composed of an adiabatically expanded composite gas 
containing a high CR component. 



In the right-hand panel of Fig. 1, we show the result for the 
case of the Mach number Al = 10 obtained with gadget-2 at time 
t = 0.3. The agreement with the analytic solution is good, while the 
discontinuities are resolved within 2-3 SPH smoothing lengths. 
Despite the shock broadening, the post-shock quantities are calcu- 
lated very accurately. In the case of composite gas, our formaUsm 
is clearly able to detect the shock and measure its strength with a 
Mach number accuracy better than 10%. Although the total pres- 
sure remains constant across the contact discontinuity, the partial 
pressure of CRs and thermal gas interestingly are changing. This 
behaviour reflects the adiabatic compression of the CR pressure 
component across the shock wave. A posteriori, this justifies our 
procedure of inferring the thermal pressure jump at the shock for a 
composite of CRs and thermal gas in equation (37). 

In the right-hand panel of Fig. 2, we show the Mach num- 
ber distributions weighted by the change of dissipated energy per 
time interval, (dMth/(d log At)) for our eight shock-tubes. While our 
formalism is able to measure the shock strength with a Mach num- 
ber accuracy better than 10%, the distributions are sharply peaked. 
This demonstrates the reliability of our formalism to measure shock 
strengths for the composite gas instantaneously during SPH simu- 
lations. 

The bottom panel shows the change of dissipated energy per 
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time interval, (!<,(,> together with the shock-injected CR energy 
(«CRjnj)- Concerning the amount of injected CR energy, we ne- 
glected cooling processes such as Coulomb interactions with ther- 
mal particles: this would effectively result in a density dependent 
recalibration of the maximum CR energy efficiency ^max of the oth- 
erwise arbitrary absolute value of our fiducial density. In the case 
of high Mach numbers, there is a good agreement between the sim- 
ulated and theoretically expected shock-injected CR energy while 
there are discrepancies at low Mach numbers: our formalism es- 
timates volume averaged Mach numbers with an accuracy better 
than 10%; this vmcertainty translates to estimates of the density 
jump and the thermal pressure jump with a scatter among 
different SPH particles. In the regime of weak shocks, the CR en- 
ergy efficiency due to diffusive shock acceleration finj is extremely 
sensitive to these two quantities, leading to larger uncertainties for 
the shock-injected CR energy in the case of a high CR load. How- 
ever, the overall trend for the shock-injected CR energy can still be 
matched in such an extreme physical environment. 



6 NON-RADIATIVE COSMOLOGICAL SIMULATIONS 
6.1 Simulation setup 

As a first application of our formalism, we are here interested in 
studying the spatial distribution of cosmological structure forma- 
tion shocks in combination with Mach number statistics. We focus 
on the "concordance" cosmological cold dark matter model with a 
cosmological constant (ACDM). The cosmological parameters of 
our model are: = Sidm + Sib = 0.3, Sib = 0.04, Qa = 0.7, 
h = 0.7, n = I, and erg = 0.9. Here, Si^ denotes the total matter 
density in units of the critical density for geometrical closure, Pait = 
3//Q/(87tG). Sib and Ha denote the densities of baryons and the cos- 
mological constant at the present day. The Hubble constant at the 
present day is parametrized as Ho = VX)h km s^'Mpc"', while n 
denotes the spectral index of the primordial power- spectrum, and 
(Tg is the rms linear mass fluctuation within a sphere of radius 
8 /!"'Mpc extrapolated to z = 0. This model yields a reasonable fit 
to current cosmological constraints and provides a good framework 
for investigating cosmological shocks. 

Our simulations were carried out with an updated and ex- 
tended version of the distributed-memory parallel TreeSPH code 
GADGET-2 (Springel 2005; Springel et al. 2001) including now 
self-consistent cosmic ray physics (EnBlin et al. 2006; Jubelgas 
et al. 2006). Our reference simulation employed 2 x 256' par- 
ticles which were simulated within a periodic box of comov- 
ing size 100 h'^Mpc, so the dark matter particles had masses of 
4.3 X 10' /j-' Mg and the SPH particles 6.6 x 10« hr^ Mq. The SPH 
densities were computed from 32 neighbours which leads to our 
minimum gas resolution of approximately 2 x lO"'/?"' Mq. The 
gravitational force softening was of a spline form (e.g., Hem- 
quist & Katz 1989) with a Plummer-equivalent softening length of 
13 /i^'kpc comoving. In order to test our numerical resolution, we 
additionally simulated the same cosmological model with 2 x 128' 
particles, with a softening length twice that of the reference simu- 
lation. 

Initial conditions were laid down by perturbing a homoge- 
neous particle distribution with a realization of a Gaussian ran- 
dom field with the ACDM linear power spectrum. The displace- 
ment field in Fourier space was constructed using the Zel'dovich 
approximation, with the amplitude of each random phase mode 
drawn from a Rayleigh distribution. For the initial redshift we chose 



1 + Zinit - 50 which translates to an initial temperature of the gas 
of Tinit = 57 K. This reflects the fact that the baryons are thermally 
coupled to the CMB photons via Compton interactions with the 
residual free electrons after the universe became transparent until it 
eventually decoupled at 1 -i- Zdec - 100(ab/jV0.0125)^''^ In all our 
simulations, we stored the full particle data at 100 output times, 
equally spaced in log(l -I- z) between z = 40 and z = 0. 

In order to investigate the effects of reionisation on the Mach 
number statistics, we additionally perform two similar simulations 
which contain a simple reionisation model where we impose a min- 
imum gas temperature of r = 10* K at a redshift of z = 10 to all 
SPH particles. We decided to adopt this simplified model to study 
its effect on the Mach number statistics rather than a more compli- 
cated reionisation history. A more realistic scenario might be to add 
energy only to gas within haloes above a certain density in combi- 
nation with energy input from QSO activity, and to describe the 
merging of the reionisation fronts and their evolution into the lower 
density regions (e.g., Ciardi et al. 2003). 

The simulation reported here follow only non-radiative gas 
physics. We neglected several physical processes, such as radia- 
tive cooling, galaxy/star formation, and feedback from galaxies and 
stars including cosmic ray pressure. Our primary focus are shocks 
that are mostly outside the cluster core regions. Thus, the conclu- 
sions drawn in this work should not be significantly weakened by 
the exclusion of these additional radiative processes. 

In contrast to the idealised shock tube experiment where all 
particles that are shocked experience the same shock strength, in 
cosmological simulations there might be a a distribution of Mach 
numbers for a given region of space because of curvature effects, 
multiple shocks, etc. Thus the Mach number estimation in shock 
tubes involves averaging over many particles all of which are expe- 
riencing a shock of a given Mach number, whereas in the cosmo- 
logical simulations the averaging has to be done over Mach number 
also. This might introduce a scatter to the Mach mmiber estimation 
of a single particle in cosmological simulations which is difficult to 
quantify. We are confident that this effect has only a minor impact 
on our results because they agree well with results of similar stud- 
ies that used Eulerian structure formation simulations (Ryu et al. 
2003). 



6.2 Visualisation of the Mach number 

In the SPH formalism, continuous fields A{x) such as the hydrody- 
namical quantities are represented by the values A, at discrete par- 
ticle positions r, = (x,, y, , z,) with a local spatial resolution given by 
the SPH smoothing length h,. To visualise a scalar quantity in two 
dimensions we employ the mass conserving scatter approach for 
the projection, where the particle's smoothing kernel is distributed 
onto cells of a Cartesian grid which is characterised by its physical 
mesh size g. The line-of-sight integration of any quantity A(x) at 
the pixel at position r = {x,y,z) is determined as the average of 
integration of all lines of sight passing through the pixel. 



i \.Jx-gl2 Jy-g/2 J -hi \»i 



,(38) 



with r = ^(xi — + (y,- - yY + z^, and where the summation is 
extended over all particles in the desired slice of projection. The 
function "K is the dimensionless spherically symmetric cubic spline 
kernel suggested by Monaghan & Lattanzio (1985). 
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Figure 3. Visualisation of a non-radiative cosmological simulation at redshift z = 2 {top panels) and z = (bottom panels). Shown are the overdensity of 
the gas {left-hand side) and the density weighed gas temperature {right-hand side). These pictures have a comovrng side length of 100 /i"' Mpc while the 
projection length along the line-of-sight amounts to 10 Mpc. 



The left-hand side of Fig. 3 shows the time evolution of the 
density contrast S averaged over the line-of-sight with a comoving 
projection length Lproj = 10 hr^ Mpc: 



(l + 5..AxS) 



los Lproj Pcnl fib 



(39) 



where E denotes the surface mass density. The fine-spun cosmic 
web at high redshift evolves into a much more pronounced, knotty 
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Figure 4. Mach number visualisation of a non-radiative cosmological simulation at redshift z = 2 (top panels) and z = (bottom panels). The colour hue of the 
maps on the left-hand side encodes the spatial Mach number distribution weighted by the rate of energy dissipation at the shocks, normalised to the simulation 
volume. The maps on the right-hand side show instead the Mach number distribution weighted by the rate of CR energy injection above q = 0.8, the threshold 
of hadronic interactions. The brightness of each pixel is determined by the respective weights, i.e. by the energy production density. These pictures have a 
comoving length of 100 /i"' Mpc on a side. Most of the energy is dissipated in weak shocks which are situated in the internal regions of groups or clusters, 
while collapsed cosmological structures are surrounded by strong external shocks (shown in blue). 
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and filamentary structure at late times, as a result of the hierarchical 
structure formation process driven by gravity. 

The right-hand side of Fig. 3 shows the time evolution of 
the density weighted temperature averaged over the line-of-sight. 
Again, the growth of galaxy clusters visible as large bright regions 
with temperatures around lO^K is clearly visible. Through dissipa- 
tion, the shock waves convert part of the gravitational energy asso- 
ciated with cosmological structure formation into internal energy of 
the gas, apart from the additional contribution due to adiabatic com- 
pression caused by the material that falls in at later times and itself 
is compressed at these shock waves. The large black regions show 
voids which cool down during cosmic evolution due to two effects: 
while the universe expands, non-relativistic gas is adiabatically ex- 
panded and cools according to T cc V'"'' cc a^^ for y = 5/3 when 
shock heating is still absent. Secondly, matter is flowing towards 
filaments during structure formation, hence the voids get depleted, 
providing an additional adiabatic expansion of the remaining mate- 
rial. 

Fig. 4 shows a visualisation of the responsible structure for- 
mation shocks and their corresponding strengths. The colour scal- 
ing represents the spatial Mach number distribution weighted by 
the rate of energy dissipation at the shocks, and normalised to the 
simulation volume (left-hand side). The Mach nvunber distribution 
weighted by the rate of CR energy injection is shown in the right- 
hand side, again normalised to the simulation volume. The bright- 
ness of these pixels scales with the respective weights, i.e. by the 
rates of energy dissipation or injection, respectively. The spatial 
Mach number distribution impressively reflects the nonlinear struc- 
tures and voids of the density and temperature maps of Fig. 3. It is 
apparent that most of the energy is dissipated in weak shocks which 
are situated in the internal regions of groups or clusters while col- 
lapsed cosmological structures are surrounded by external strong 
shocks (shown in blue). These external shocks are often referred 
to as 'first shocks', because here the compressed gas has been pro- 
cessed for the first time in its cosmic history through shock waves. 

Following Ryu et al. (2003), we classify structure formation 
shocks into two broad populations which are labelled as internal 
and external shocks, depending on whether or not the associated 
pre-shock gas was previously shocked. Rather than using a ther- 
modynamical criterion such as the temperature, we prefer a crite- 
rion such as the overdensity 6 in order not to confuse the shock 
definition once we will follow radiatively cooling gas in galaxies 
(in practice, we use the criterion of a critical pre-shock overdensity 
6 > \0 for the classification of an internal shock). External shocks 
surround filaments, sheets, and haloes while internal shocks are lo- 
cated within the regions bound by external shocks and are created 
by flow motions accompanying hierarchical structure formation. 
For more detailed studies, internal shocks can be further divided 
into three types of shock waves: (1) accretion shocks caused by in- 
falling gas from sheets to filaments or haloes and from filaments to 
haloes, (2) merger shocks resulting from merging haloes, and (3) 
supersonic chaotic flow shocks inside nonlinear structures which 
are produced in the course of hierarchical clustering. 

In contrast to the present time, the comoving surface area of 
external shock waves surpasses that of internal shocks at high red- 
shift, due to the small fraction of mass bound in large haloes and 
the simultaneous existence of an all pervading fine-spun cosmic 
web with large surface area. Also, there the thermal gas has a low 
sound velocity c - ^fyPjp - ^jyiy - 1)k owing to the low temper- 
ature, so once the dilfuse gas breaks on mildly nonlinear structures, 
strong shock waves develop that are characterised by high Mach 
numbers M = f s/c. Nevertheless, the energy dissipation rate in in- 



Mach number M 




Figure 5. CR energy injection efficiency for the diffusive shock acceleration 
process. Shown is the CR energy injection efficiency fi„j = ecR/£diss tor 
the ffiree post-shock temperatures kT2/key = 0.01,0.3, and 10. We inject 
only CRs above the kinematic threshold ^threshold = 0.83 of the hadronic 
CRp-p interaction which are able to produce pions that successively decay 
into secondary electrons, neutrinos, and y-rays. We choose equipartition 
between injected CR energy and dissipated thermal non-relativistic energy 
as saturation value of the CR energy injection efficiency, i.e. fmax = 0.5 
(Ryu et al. 2003). 

temal shocks is always higher compared to external shocks because 
the mean shock speed and pre-shock gas densities are significantly 
larger for internal shocks. 

We use the same colour and brightness scale for the Mach 
number distribution weighted by the injected CR energy rate nor- 
maUsed to the simulation volume (right-hand side of Fig. 4). We 
emphasize two important points which have fundamental implica- 
tions for the CR population in galaxy clusters. (1) There is an ab- 
sence of weak shocks (shown in yellow) when the Mach number 
distribution is weighted by the injected CR energy. This reflects the 
Mach number dependent energy injection efficiency: the CR injec- 
tion is saturated for strong shocks which leads to similar spatial dis- 
tribution of both weightings, by dissipated energy as well as by in- 
jected CR energy. In contrast, most of the dissipated energy is ther- 
malized in weak shocks and only small parts are used for the accel- 
eration of relativistic particles (compare Fig. 5). (2) The mechanism 
of energy dissipation at shocks is very density dependent, implying 
a tight correlation of weak internal shocks and the amount of dis- 
sipated energy. This can be seen by the strongly peaked brightness 
distribution of the dissipated energy rate towards the cluster cen- 
tres. For the CR-weighted case, this correlation is counteracted by 
the CR energy injection efficiency leading to a smoother brightness 
distribution of the CR energy injection. This has the important im- 
plication that the ratio of CR injected energy to dissipated thermal 
energy is increasing as the density declines. Relative to the thermal 
non-relativistic energy density, the CR energy density is dynami- 
cally more important at the outer cluster regions and dynamically 
less important at the cluster centres. 

6.3 Mach number statistics 

6.3.1 Influence of reionisation 

To quantify previous considerations, we compute the differential 
Mach mmiber distribution weighted by the dissipated energy nor- 
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malised to the simulation volume d^e(iiss(a, M)/(d log a d log M) at 
different redshifts. The top left-hand panel of Fig. 6 shows this 
Mach number distribution for our reference simulation with reion- 
isation (showing a resolution of 2 x 256'), while the top right- 
hand panel shows this distribution for the simulation without reion- 
isation. The lower left-hand panel shows both distributions inte- 
grated over the scale factor, dejiss(AI)/ (d log M), in addition to the 
Mach number distribution weighted by the injected CR energy nor- 
malised to the simulation volume, d£cR(Wt)/(dlog Al) (see Sec- 
tion 6.3.2). Internal shocks are shown with dotted lines and exter- 
nal shocks with dashed lines. The lower right-hand panel shows the 
evolution of the dissipated energy due to shock waves with scale 
factor, d£diss(a)/(dloga), for the models with and without reionisa- 
tion. 

Several important points are apparent. (1) The median of the 
Mach munber distribution weighted by the dissipated energy de- 
creases as cosmic time evolves, i.e. the average shock becomes 
weaker at later times. (2) There is an increasing amount of en- 
ergy dissipated at shock waves as the universe evolves because the 
mean shock speed is significantly growing when the characteris- 
tic mass becomes larger with time. This trend starts to level off at 
redshift z - I although the median Mach number in shocks con- 
tinuous to decrease. (3) Reionisation influences the Mach number 
distribution predominantly at early times (however after reionisa- 
tion took place) and suppresses strong external shock waves effi- 
ciently. The reason is that reionisation of the thermal gas increases 
its sound speed c = -\JynkT/p dramatically, so that weaker shocks 
are produced for the same shock velocities. (4) The time integrated 
Mach number distribution weighted by the dissipated energy peaks 
at Mach numbers 1 < A1 < 3. The main contribution in terms 
of energy dissipation originates from internal shocks because of 
enhanced pre-shock densities and mean shock speeds. (5) Exter- 
nal shocks dominate the Mach number distribution at early times 
while internal shocks take over at z ^ 9 (depending somewhat on 
the resolution of the simulation). Their amount of dissipated energy 
surpasses that in external shocks by over an order of magnitude at 
the present time. Internal shocks play a more important role than 
external shocks in dissipating energy associated with structure for- 
mation. 

The total shock-dissipated energy in our simulation box 
amounts to E^i^s = 2.27 x 10*"* erg. This translates to a mean en- 
ergy deposition per particle of £diss /"/(Pcrit^^b ^) = 0.63 keV, where 
IX = 4OTp/(3 + 5Xii) is the mean particle weight assuming full ioni- 
sation and = 0.76 is the primordial hydrogen mass fraction. Our 
results agree well with those of Ryu et al. (2003) in the case of in- 
ternal shocks while our external shocks tend to be weaker. This can 
be attributed to our differing definition of internal/external shocks 
as we prefer a density criterion and use the critical pre-shock over- 
density 6 > 10 for the classification of an internal shock. 

6.3.2 Cosmic ray acceleration 

In our non-radiative cosmological simulations we additionally cal- 
culate the expected shock-injected CR energy using our formalism 
of diffusive shock acceleration described in EnlSlin et al. (2006). 
This formalism follows a model based on plasma physics for 
the leakage of thermal ions into the CR population. However, 
in the present analysis, the injected CR energy dMcR.inj /(d log a) 
was only monitored and not dynamically tracked. In our model, 
the CR population is described by single power-law distribution 
which is uniquely determined by the dimensionless momentum 
cutoff q, the normalisation C, and the spectral index a. Consid- 



ering shock injected CRs only, the spectral index is determined by 
a = {Xs + 2)/{Xs - 1), where Xs denotes the density jump at the 
shock. 

Our simplified model for the diffusive shock acceleration fails 
in the limit of weak shocks and over-predicts the injection effi- 
ciency. Especially in this regime. Coulomb losses have to be taken 
into account which remove the low-energetic part of the injected 
CR spectrum efficiently on a short timescale giving rise to an ef- 
fective CR energy efficiency. Thus, the instantaneous injected CR 
energy dMcR,inj /(d log a) depends on the simulation timestep and the 
resolution. To provide a resolution independent statement about the 
injected CR energy, we decided to rethermalise the injected CR 
energy below the cutoff ^threshold = 0.83. This cutoff has the de- 
sired property, that it coincides with the kinematic threshold of the 
hadronic CR p-p interaction to produce pions which decay into sec- 
ondary electrons (and neutrinos) and y-rays: 

JU* + Vjj/Vij ^ -I- Ve/Ve + Vf, + Vjj 

2y. 

Only CR protons above this kinematic threshold are therefore visi- 
ble through their decay products in both the y-ray and radio bands 
via radiative processes, making them directly observationally de- 
tectable. 

The lower left-hand panel of Fig. 6 shows the Mach nvrni- 
ber distribution weighted by the injected CR energy rate and 
normalised to the simulation volume, decR(A1)/(dlog Al) (solid 
green). The effect of the CR injection efficiency = ecR/^diss 
can easily be seen: while the CR injection is saturated for strong 
shocks to (fmax = 0.5, in weak shocks most of the dissipated en- 
ergy is thermalized and only small parts are used for the acceler- 
ation of relativistic particles. Effectively, this shifts the maximum 
and the mean value of the Mach number distribution weighted by 
the shock-dissipated energy towards higher values in the case of 
the distribution weighted by the injected CR energy. This effect 
is even stronger when considering only CRs with a lower cut- 
off" ^ = 10, 30 which are responsible for radio haloes observed 
at frequencies above 100 MHz, assuming typical magnetic field 
strengths of S = 10, 1 jjG, respectively. This follows from the 
mono-energetic approximation of the hadronic electron production 
and synchrotron formula, 

3eB , q m„ 

n = 7^ ri, where 75-77— (40) 

InnieC 16 fWe 

and e denotes the elementary charge. 

As the regime of strong shocks is dominated by external 
shocks where the CR injection is saturated, CRs are dynamically 
more important in dilute regions and dynamically less important 
at the cluster centres compared to the thermal non-relativistic gas. 
As weak shocks are mainly internal shocks we have to distinguish 
between their different appearance: strong internal shocks are most 
probably accretion shocks produced by infalling gas from sheets 
or filaments towards clusters, or peripheral merger shocks which 
steepen as they propagate outwards in the shallow cluster potential, 
highlighting the importance of CR injection in the outer cluster re- 
gions relative to thermally dissipated gas at shocks. In contrast, CR 
injection is dynanucally less important in the case of flow shocks 
at the cluster centres or merging shock waves traversing the clus- 
ter centre. From these considerations we again draw the important 
conclusion that the ratio of CR injected energy to dissipated ther- 
mal energy at shocks is an increasing function of decreasing den- 
sity. Such a CR distribution is required within galaxy clusters to ex- 
plain the diffuse radio synchrotron emission of galaxy clusters (so- 
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Figure 6. Influence of reionisation on the Mach number statistics of non-radiative cosmological simulations. The top left-hand panel shows the differential 
Mach number distribution d^S(jiss(a, A1)/(d log a d log M) for our reference simulation with reionisation while the top right-hand panel shows this distribution 
for the simulation without reionisation. The lower left-hand panel shows both distributions integrated over the scale factor, de(iiss(A1)/(d log M) in addition to 
the Mach number distribution weighted by the injected CR energy rate normalised to the simulation volume, d£cR(A1)/(dlog M) (green). Internal shocks are 
shown with dotted lines and external shocks with dashed lines. The lower right-hand panel shows the evolution of the dissipated energy due to shock waves 
with scale factor, d£(iiss(«)/(dloga). The models with and without reionisation lie on top of each other. 



called radio haloes) within the hadronic model of secondary elec- 
trons. For that, we assume a stationary CR electron spectrum which 
balances hadronic injection of secondaries and synchrotron and 
inverse Compton cooling processes (Brunetti 2002; Pfrommer & 
EnBlin 2004a,b). However, to make more precise statements about 
the origin of cluster radio haloes, more work is needed which stud- 
ies the effect of the CR dynamics including CR diffusion and other 
CR injection processes such as supernovae driven galactic winds. 



6..?..? Resolution study 

To study numerical convergence we perform two additional simu- 
lations with a resolution of 2 x 128^^, respectively, for our models 
with and without reionisation. Fig. 7 shows this resolution study for 
non-radiative cosmological simulations with a reionisation epoch at 
z = 10. The lower right-hand panel of Fig. 7 shows the evolution 
of the density-weighted temperature with redshift, <(1 -I- S)T)v (z). 
Shown are different resolutions in our models with and without 
reionisation. The two differently resolved simulations (2x256'' and 
2 X 128-') have converged well at redshifts z < 4. In our reference 
simulation, the adiabatic decay of the mean temperature is halted at 
slightly higher redshift: because of the better mass resolution of this 
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Figure 7. Resolution study: Mach number statistics for non-radiative cosmological simulations with a reionisation epoch at z = 10. The top left-hand panel 
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evolution of the density-weighted temperature with redshift. Shown are different resolutions in our models with and without reionisation. 



simulation, nonlinear structures of smaller mass can be resolved 
earlier while converting part of their gravitational binding energy 
into internal energy through structure formation shock waves. In 
the simulation with reionisation, the temperature increases discon- 
tinuously at z = 10, declines again with the adiabatic expansion, 
until shock heating takes over at z ~ 7-8 (depending on the 
resolution of the simulation). At z = 0, all simulations yield a 
mean density-weighted temperature of ^ 0.3 keV. Comparing this 
density-weighted energy to the shock-deposited mean energy per 
particle of ^diss - 0.63 keV, we obtain the mean adiabatic com- 
pression factor of the cosmic plasma, {A;r/[(7-l)£'(]iss]j''<''"'' ^ 0.6. 
After the plasma has been shock-heated, relaxation processes in the 
course of virialisation let the plasma expand adiabatically on aver- 



age. Secondly, mildly non-linear structures characterised by a shal- 
low gravitational potential are partly effected by the Hubble flow 
which forces them to adiabatically expand. 

The top left-hand panel shows the differential Mach number 
distribution d^ediss(a, A1)/(dlogadlog AI) for our reference sim- 
ulation with a resolution of 2 x 256' while the top right-hand 
panel shows this distribution for the simulation with a resolution 
of 2 X 128'. The lower left-hand panel shows both distributions in- 
tegrated over the scale factor, d£:diss(A1)/(dlog AI). Internal shocks 
are shown with dotted lines and external shocks with dashed lines. 
One immediately realizes that the question if the first shocks are 
fully converged among simulations of different resolution is not 
well posed because nonlinear structures of smaller mass can be re- 
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solved collapsing earlier in higher resolution simulations. Accord- 
ingly, the differential Mach number distribution is not well con- 
verged at redshifts z > 6 while the distribution is well converged 
for z < 3. Since most of the energy is dissipated at late times, where 
our differential Mach number distribution is well converged, the in- 
tegrated distribution dediss(Al)/(dlog Al) shows only marginal dif- 
ferences among the differently resolved simulations. In particular, 
our statements about CR injection at structure formation shocks are 
robust with respect to resolution issues. 



7 SUMMARY AND CONCLUSIONS 

We provide a formalism for identifying and estimating the strength 
of structure formation shocks in cosmological SPH simulations on- 
the-fly, both for non-relativistic thermal gas as well as for a plasma 
composed of a mixture of cosmic rays (CRs) and thermal gas. In 
addition, we derive an analytical solution to the one-dimensional 
Riemann shock tube problem for the composite plasma of CRs and 
thermal gas (Appendix B). In the case of non-relativistic thermal 
gas, shock-tube simulations within a periodic three-dimensional 
box that is longer in x-direction than in the other two dimensions 
show that our formalism is able to unambiguously detect and accu- 
rately measure the Mach numbers of shocks, while in the case of 
plasma composed of cosmic rays (CRs) and thermal gas, the Mach 
numbers of shocks are estimated with an accuracy better than 10%. 
In both cases, we find a very good agreement of the averaged simu- 
lated hydrodynamical quantities (such as density, pressure, and ve- 
locity) and the analytical solutions. Using our formalism for diffu- 
sive shock acceleration, we additionally calculate and monitor the 
shock-injected CR energy, but without dynamically tracking this 
CR energy component; the latter will be studied in forthcoming 
work. The good agreement between the simulated and theoretically 
expected shock-injected CR energy demonstrates that our formal- 
ism is reliably able to accelerate CRs instantaneously during the 
simulation. 

Subsequently, we identified and studied structure formation 
shock waves using cosmological N-body/hydrodynamical SPH 
simulations for a concordance ACDM universe in a periodic cu- 
bic box of comoving size lOOA^'Mpc. We performed simulations 
with and without a reionisation epoch at z = 10 in order to inves- 
tigate the effects of reionisation on the Mach number distribution. 
Our sets of simulations follow only non-radiative gas physics where 
we neglected additional physical processes, such as radiative cool- 
ing, star formation, and feedback from galaxies and stars includ- 
ing cosmic ray pressure. Since we are mainly interested in shock 
waves situated mostly outside the cluster core regions, the conclu- 
sions drawn in this article should not be significantly weakened 
by the exclusion of those radiative processes. Furthermore, these 
simplifications align our work with the mesh-based simulations of 
Ryu et al. (2003) and enable a direct comparison and verification 
of our results. We classify cosmological shock waves as internal 
and external shocks, depending on whether or not the associated 
pre-shock gas was previously shocked (cf. Ryu et al. 2003). Rather 
than using a thermodynamical criterion such as the temperature, we 
prefer a density criterion such as the overdensity 6 in order not to 
confuse the shock definition once we will follow radiatively cool- 
ing gas in galaxies. External shocks surround filaments, sheets, and 
haloes where the pristine adiabatically cooling gas is shocked for 
the first time. Internal shocks on the other hand are located within 
the regions bound by external shocks and are created by flow mo- 
tions accompanying hierarchical structure formation. Their popula- 



tion includes accretion shocks produced by infalling material along 
the filaments into clusters, merger shocks resulting from infalling 
haloes, and flow shocks inside nonUnear structures which are ex- 
cited by supersonic motions of subclumps. 

The Mach number distribution weighted by the dissipated en- 
ergy shows in detail that most of the energy is dissipated in weak 
shocks which are situated in the internal regions of groups or clus- 
ters while collapsed cosmological structures are surrounded by ex- 
ternal strong shocks which have a minor impact on the energy bal- 
ance. The evolution of the Mach number distribution shows that 
the average shock strength becomes weaker at later times while 
there is an increasing amount of energy dissipated at shock waves 
as cosmic time evolves because the mean shock speed increases 
together with the characteristic mass of haloes forming during cos- 
mic structure formation. For the same reason, internal shocks play 
a more important role than external shocks in dissipating energy as- 
sociated with structure formation, especially at small redshift. The 
energy input through reionisation processes influences the Mach 
number distribution primarily during a period following the reion- 
isation era and suppresses strong external shock waves efficiently 
because of the significant increase of the sound speed of the inter- 
galactic medium. 

Weighting the Mach number distribution by the injected CR 
energy shows the potential dynamical implications of CR popula- 
tions in galaxy clusters and haloes: the maximum and the mean 
value of the Mach number distribution, weighted by the shock- 
dissipated energy, is effectively shifted towards higher values of the 
Mach number when the distribution is weighted by the injected CR 
energy. In other words, the average shock wave responsible for CR 
energy injection is stronger compared to the average shock which 
thermalizes the plasma. The fundamental reason for this lies in the 
theory of diffusive shock acceleration at coUisionless shock waves 
and can be phenomenologically described by a CR injection effi- 
ciency: while the CR injection is saturated to an almost equiparti- 
tion value between injected CR energy and dissipated thermal en- 
ergy for strong shocks, in weak shocks most of the dissipated en- 
ergy is thermalized and only small parts are used for the accelera- 
tion of relativistic particles. Relative to the thermal non-relativistic 
energy density, the shock-injected CR energy density is dynami- 
cally more important at the outer dilute cluster regions and less 
important at the cluster centres since weak shock waves predom- 
inantly occur in high-density regions. This has the crucial conse- 
quence that the ratio of CR injected energy to dissipated thermal 
energy is an increasing function as the density decUnes. Such a 
CR distribution within galaxy cluster is required to explain the dif- 
fuse radio synchrotron emission of galaxy clusters (so-called radio 
haloes) within the hadronic model of secondary electrons. In or- 
der to draw thorough conclusions about the origin of cluster radio 
haloes, more work is needed which studies the effect of the CR dy- 
namics comprising of CR injection and cooling processes as well 
as CR diffusion mechanisms. 

We note that our new formalism for shock-detection in SPH 
simulations should have a range of interesting applications in sim- 
ulations of galaxy formation. For example, when combined with 
radiative dissipation and star formation, our method can be used 
to study CR injection by supernova shocks, or to construct models 
for shock-induced star formation in the interstellar medium (e.g. 
Barnes 2004). It should also be useful to improve the accuracy of 
predictions for the production of -y-rays by intergalactic shocks (e.g. 
Keshet et al. 2003). 
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APPENDIX A: RIEMANN SHOCK-TUBE PROBLEM 



The Riemann shock-tube calculation of Sod (1978) has become a generally accepted test of numerical hydrodynamical codes. As a baseline 
for later extension, we present in the section the quasi-analytical solution for the Riemann problem in the standard case of a polytropic 
gas. Then, in Appendix B we derive the quasi-analytic solution in the case of a gas composed of CRs and thermal gas, where the effective 
adiabatic index depends on the different equations of state and changes across the shock-tube. 

In the following, we summarise the key considerations which lead to the solution of the Riemann problem, for completeness (see e.g. 
Courant & Friedrichs 1948; Toro 1997; Rasio & Shapiro 1991, for a compact representation). For the initial state, we assume a state with 
higher pressure in the left-half space without loss of generaUty. At any time r > 0, this leads to the development of five regions of gas with 
different hydrodynamical states which are numbered in ascending order from the right-hand side. These regions are separated by the head 
and the tail of the leftwards propagating rarefaction wave, and the rightwards propagating contact discontinuity and the shock wave. Mass, 
momentum and energy conservation laws are represented by the generalised Rankine-Hugoniot conditions for a given coordinate system: 



Va\p\ 



= Ipv], 

= Ipv^+P], 



(Al) 



Py +e + P|w 



Here Vi denotes the speed of the discontinuity under consideration with respect to our coordinate system and we introduced the abbreviation 
[F] = F; - Fj for the jump of some quantity F across the discontinuity. Within the leftwards propagating rarefaction wave, the generalised 
Riemann invariants yield an isentropic change of state, d* = 0, and conserve the quantity F*: 



Jo P' 



2c(p) 

r-1 



: const. 



(A2) 



For the last step, we assumed a polytropic equation of state P = Ap^ . Appropriately combining these equations, the solution can be expressed 
as follows: 



p(x, = 



P5, 

P5[V 
P3. 
P2, 
Pi. 



(I-P^)] 



2/(r-i) 



v(x, t) = 



^5, 



|2r/(r-i) 



Pi 
Pu 



-c^t < X<, -ft?, 
-Vtt < X<, V2t, 
V2t < X < Vst, 
X > Vst, 

X < -Cst, 
-Cst < X < -Vtt, 
-Vtt < X <: Vst, 
X > Vst, 



(A3) 



0, X < -cst, 

(1 -/i^)(j +C5), -C5t<X<-Vtt, 

V2 = V3, -Vtt < X < Vst, 

0, X> Vst. 



(A4) 



(A5) 



Here yu^ = (7 - l)/(y -1- 1), Ci = ^jyPilpi, and C5 = ^fyPTfpl are the speeds of sound, Vi is the speed of propagation of the rarefaction wave's 
tail, and Vs is the shock speed. The post-shock pressure is obtained by solving (numerically) the non-linear equation, which is derived from 
the Rankine-Hugoniot conditions over the shock while ensuring the conservation of the two Riemann invariants of equation (A2): 



1 



C5 



y{P2lPi+ff) (y-l)ci 



1 - 



(r-i)/(2r) 



0. 



(A6) 



The density on the left-hand side of the contact discontinuity is ps = psiPilPs)^^'' , since the gas is adiabatically connected to the left-hand 
side. The post-shock density p2 is also derived from the Rankine-Hugoniot conditions. 



P2 = Pi 



(P2+p'Pi 



{Pl+H^P2} 

The post-shock gas velocity V2 is obtained from the rarefaction wave equation, xjt = v — c, and usage of the Riemann invariant F*: 

2C5 



V2=V3 



(r-1) 



p \(r-i)/(2r) 



(A7) 



(AS) 



and from equation (A5) we derive the speed of propagation of the rarefaction wave's tail i/t = C5 - 1/2/(1 - j"^)- Finally, mass conservation 
across the shock yields 

V2 



1 -P1/P2' 



(A9) 
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APPENDIX B: RIEMANN SHOCK-TUBE PROBLEM FOR A COMPOSITE OF COSMIC RAYS AND THERMAL GAS 
Bl Derivation 

In contrast to the previous case, the composite of CRs and thermal gas does not obey a polytropic equation of state. In this section, we 
present an analytical derivation of the Riemann shock-tube problem for the composite of polytropic gas and a component that is adiabatically 
compressed at the shock such as relativistic gas or a homogeneous magnetic field which is parallel to the shock front. For the analytical 
derivation, we adopt the following two approximations: (i) we assume the CR adiabatic index (equation 6) to be constant over the shock- 
tube, and (ii) we neglect CR diffusion. The first assumption is justified as long as the CR pressure is not dominated by trans-relativistic CRs 
of low energy while the second assumption is a strong simplification with respect to simulating realistic shocks including CRs (Kang & Jones 
2005). However, including CR diffusion complicates the problem significantly such that it is not any more analytically tractable. 

For the initial state, we again assume a state with higher pressure in the left-half space. At any time f > 0, five regions of gas with 
different hydrodynamical states coexist, and are numbered in ascending order from the right-hand side. We use the notation P, = Pcr.i + Ah,/ 
and Sj = EcR,/ + ^thj for the composite quantities in region /. The full solution of the initial value problem consists of determining 12 unknown 
quantities in the regions (2) and (3): p2, V2, Pcri, Pthi, Scr2, sm, and p3, vj, Pqrs, Pm, £cr3, Sm- The thermal gas obeys a polytropic equation 
of state, i.e. eth,; = ^th,i7(7th - 1) for i € (2, 3) and the regions (2) and (3) are separated by a contact discontinuity, implying vanishing mass 
flux through it and thus, V2 = 1/3 and P2 = P3. This reduces the dimensionality of our problem to 8 imknowns. In our approximation, the CRs 
are adiabatically expanded over the rarefaction wave and adiabatically compressed at the shock while obeying a pol5?tropic equation of state: 



'(pf) ' = ecRi(^) 



rcR 



(Bl) 



which further reduces the dimensionality by 4 unknowns. Moreover, the thermal gas is also adiabatically expanded over the rarefaction wave 
yielding Pth3 = PihiiPil Piy^ ■ Hence, we need 3 more linearly independent equations for the solution: 2 are obtained by considering the 
Rankine-Hugoniot conditions (equation Al) in a stationary system of reference with Vi = Vg. The last equation is given by the Riemann 
invariant F*, where the effective speed of sound is given by c = ^fy^gP/P' 



^^^dp' = u + Hp) = const, with I(p) = f Jacrx''cr-3 + Athxy^^-^dx. (B2) 
P' Jo ^ 



Here, we use the abbreviations A,- = yjA,- where i € (th, CR) and A- = P-p ^' denotes the invariant adiabatic function over the rarefaction 
wave. Introducing the difference of the adiabatic indices of the two populations. Ay = jtb - 7cr, the solution to the integral /(p) is given by 



I(p) = ^r^\ s with x(p} = , , . (B3) 

^7 \ Aft, I 2Ay 2Ay j AcrP^cr h- Aftprr 



sJl^'-^\ with x(o)=^^ 



nrth 



Although the second argument of the incomplete Beta-function is always negative, /(p) is well defined as long as we consider a non-zero CR 
pressure which is characterised by Aqr > 0, and ycR sufficiently far from yth, i.e. Ay > 0. For Acr = 0, the integral can be solved in closed 
form, yielding I(p) = 2c(p)/(ya, - 1). 



B2 Solution of the Riemann problem 

The densities leftwards and rightwards of the contact discontinuity, pa and p2, are obtained by solving (nvmierically) the following non-linear 
system of equations. It is derived from matching the possible post-shock states (pressure and density) with the possible post-rarefaction wave 
states while simultaneously ensuring the conservation laws over the rarefaction wave and the shock: 

MXs,Xr) = [P2(x,)-P,](x,-l)-ptx,[I(ps)-l(x,ps)f = 0, 

/2(X„ X,) = [P2(X,) + Pi] (X, - 1) + 2U,£, - S2(X,,X,)] = 0. 

Here we introduced the shock compression ratio x^ = P2/P1 and the rarefaction ratio = P3/P5. Furthermore, the implicit dependences on 
Xs and Xt can explicitly be expressed as follows, 

P2(X.) = P3{Xr) = PcR5X^''^+Pmxl'', (B5) 
PcR2{Xs) = i'cRlX^^ (B6) 

S2(x„Xr) = SCRIXI"^ + —^[P2{Xr} - PcR2iXs)l (B7) 

rth - 1 

The roots of the non-linear system of equations (equation B4) immediately yield the post-shock pressure of the fluid via equation (B5). The 
post-shock velocity 1/2 = f 3 and the shock speed Vg are then obtained from the Rankine-Hugoniot relations. 



V2 = J[P2iXr)-Pl]^^-^, (B8) 
P2P1 

(B9) 

P2 -Pl 
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Using the previous results, we can construct the solution to the generalised Riemann problem for CRs and thermal gas as follows: 



p{x, t) - 



P(x, t) = 



v(x, f) ■■ 



Ps, X < -Cjt, 

p(x, t), -C^t < X < -Vit, 

P3, -V{t < X<, V2t, 

P2, V2t < X < Vgt, 

pi, X > Vgl, 

P5, X < -c^t, 

AcR p(x, + Aflip(x, tyo-, -est <x< -Vtt, 

P2 = P3, -Vtt < X < Vgt, 

Pi, X > Vst, 

0, X < -csf, 

7 + yjAcRpix, tycR-i + Aa,p{x, 0™"^ -c^r < x< -v^t, 

V2 = V3, -Vtt < X ^ Vst, 

0, X> Vst. 



(BIO) 



(Bll) 



(B12) 



Here C5 = ■sJyeesPs/Ps is the effective speed of sound, Vt is the speed of propagation of the rarefaction wave's tail, and Vs is the shock speed. 
Matching the rarefaction wave equation to the density of the post-contact discontinuity yields Vf. 



Vt = I(P3) - /(Ps) + a/WT^^'+Z^pF'- 



(B13) 



The density within the rarefaction regime is obtained by solving (numerically) the non-linear equation for a given (x, t), which is derived 
from the rarefaction wave equation. 



I\p{x, f)] - /(ps) + y + ^Ac!Lp{x, t)y^-^ + Athp{x, 0^-1 = 0. 



(B14) 
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